#How Do Observers Assess Resolve?
#July 7, 2018
#BJPS Replication 2.R

#This file contains the replication code needed to replicate the conjoint analysis from the main text (located in the "Results and Discussion" section)

### To generate the necessary dataframes, either run BJPS Replication 1.R, or load the following saved object:
library(here) #Avoids needing to setwd()
dat3 <- get(load(file="Resolve Conjoint Cleaned No US 070718.RData"))

###### Main analysis ####

#Information about experimental results:
nrow(dat3)/2 #8090 choice tasks
length(unique(dat3$ID_numeric)) #1995 participants
nrow(dat3) #16180 country profiles

#Do parallel processing - speeds things up
library(snow)
cl <- makeCluster(8,"SOCK")
set.seed(43215)

clusterBootS2 <- function(data){
	i <- sample(unique(data$ID_numeric),length(unique(data$ID_numeric)),replace=TRUE)
	row.nums <- NULL
	for (j in 1:length(i)){
		row.nums <- c(row.nums, which(data$ID_numeric==i[j]))
	}
	return(data[row.nums,])
}

######## Main model ########

#Main analysis 
bootConjoint <- function(...){
	temp <- clusterBootS2(dat3)
	mod.temp <- lm(Chosen ~  Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, data=temp)
	return(coef(mod.temp))
}

#Weighted analysis
bootConjointW <- function(...){
	temp <- clusterBootS2(dat3)
	mod.temp <- lm(Chosen ~  Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, weights=eWeight, data=temp)
	return(coef(mod.temp))
}

clusterExport(cl, list("dat3", "bootConjoint", "bootConjointW", "clusterBootS2"))

#First: main analysis
system.time(boot.full2 <- parSapply(cl, rep(1,1500), bootConjoint)) #Main results

#Extract quantities of interest
plot.mat2 <- cbind(apply(boot.full2, 1, mean), apply(boot.full2, 1, quantile, c(0.025)), apply(boot.full2, 1, quantile, c(0.05)), apply(boot.full2, 1, quantile, c(0.95)), apply(boot.full2,1,quantile, c(0.975)))[-1,]

#Now add baseline categories
plot.mat2b <- rbind(rep(0,5), plot.mat2[1,], rep(0,5), plot.mat2[2,], rep(0,5), plot.mat2[3:4,], rep(0,5), plot.mat2[5,], rep(0,5), plot.mat2[6,], rep(0,5), plot.mat2[c(8,7),], rep(0,5), plot.mat2[9,], rep(0,5), plot.mat2[10,], rep(0,5), plot.mat2[11,],rep(0,5), plot.mat2[c(14,13,12),], rep(0,5),plot.mat2[15:16,])

#Now: weighted analysis
system.time(boot.fullW2 <- parSapply(cl, rep(1,1500), bootConjointW)) #Weighted results

#Extract quantities of interest
plot.matW <- cbind(apply(boot.fullW2, 1, mean), apply(boot.fullW2, 1, quantile, c(0.025)), apply(boot.fullW2, 1, quantile, c(0.05)), apply(boot.fullW2, 1, quantile, c(0.95)), apply(boot.fullW2,1,quantile, c(0.975)))[-1,]

#Now add baseline categories
plot.matWb <- rbind(rep(NA,5), plot.matW[1,], rep(NA,5), plot.matW[2,], rep(NA,5), plot.matW[3:4,], rep(NA,5), plot.matW[5,], rep(NA,5), plot.matW[6,], rep(NA,5), plot.matW[c(8,7),], rep(NA,5), plot.matW[9,], rep(NA,5), plot.matW[10,], rep(NA,5), plot.matW[11,],rep(NA,5), plot.matW[c(14,13,12),], rep(NA,5),plot.matW[15:16,])

#Labels for plot
textLab <- c("Low Capabilities", "High Capabilities", "Low Stakes", "High Stakes", "Dictatorship", "Democracy", "Mixed", "Ally", "Adversary", "Established Leader", "New Leader", "No Military Service", "Some Military Service", "Long Military Service", "Female Leader", "Male Leader",  "Target", "Initiator", "Against ally", "Against adversary", "Different leader, stood firm", "Same leader, stood firm", "Same leader, backed down", "Different leader, backed down", "Nothing", "Mobilized troops", "Public threat")

##### Figure 2 ######
#Clean up labels, add legend in post-production
dev.new(height=9, width=6.5)
par(mar=c(5.1,1.1,4.1,0), oma=c(0,4,0,0))
plot(plot.mat2b[,1], nrow(plot.mat2b):1, type="n", axes=FALSE, xlab="Average marginal component effect (AMCE)", ylab="",xlim=c(-0.3,0.2))

for (i in 1:nrow(plot.mat2b)){
	points(plot.mat2b[i,1], nrow(plot.mat2b)-(i-1), pch=16)
	points(plot.matWb[i,1], nrow(plot.matWb)-(i-0.75), pch=12)
	lines(plot.mat2b[i,c(2,5)], rep(nrow(plot.mat2b)-(i-1),2))
	lines(plot.matWb[i,c(2,5)], rep(nrow(plot.matWb)-(i-0.75),2))
	lines(plot.mat2b[i,c(3,4)], rep(nrow(plot.mat2b)-(i-1),2), lwd=2)
	lines(plot.matWb[i,c(3,4)], rep(nrow(plot.matWb)-(i-0.75),2), lwd=2)					
	text(-0.3, nrow(plot.mat2b)-(i-1), textLab[i], cex=0.7)			
	}
abline(v=0)
axis(1)

#Quantities of interest discussed in text:
round(cbind(plot.mat2[,1], plot.matW[,1]),digits=3)[1,] #Capabilities: +16.4% (+16.1% weighted)
round(cbind(plot.mat2[,1], plot.matW[,1]),digits=3)[2,] #Stakes: +14.8% (+14.3% weighted)
round(cbind(plot.mat2[,1], plot.matW[,1]),digits=3)[3,] #Democracy: -4.2% (-3.3% weighted)
round(cbind(plot.mat2[,1], plot.matW[,1]),digits=3)[4,] #Mixed regime: -1.4% (-1.6% weighted)
round(cbind(plot.mat2[,1], plot.matW[,1]),digits=3)[6,] #New leader: -1.5% (-2.4% weighted)
round(cbind(plot.mat2[,1], plot.matW[,1]),digits=3)[9,] #Male leader: +0.04% (+1.9% weighted)
round(cbind(plot.mat2[,1], plot.matW[,1]),digits=3)[7,] #Long military service: +7.8% (+7.9% weighted)
round(cbind(plot.mat2[,1], plot.matW[,1]),digits=3)[8,] #Some military service: +3.3% (+3.0% weighted)
round(cbind(plot.mat2[,1], plot.matW[,1]),digits=3)[13,] #Same leader, backed down: -11.1% (-13.6% weighted)
round(cbind(plot.mat2[,1], plot.matW[,1]),digits=3)[12,] #Different leader, backed down: -9.1% (-9.5% weighted)
round(cbind(plot.mat2[,1], plot.matW[,1]),digits=3)[14,] #Same leader, stood firm: +4.9% (+5.2% weighted)

######### Interactive models #######

#### Current calculus hypothesis
press1 <- dat3[,-c(2:5,8:13,15:21,23:25)] #Slimmer dataframe

bootPress1 <- function(...){
	temp <- clusterBootS2(press1)
	mod.temp <- lm(Chosen ~ backCapInt, data=temp)
	return(coef(mod.temp))
}

bootPress1W <- function(...){
	temp <- clusterBootS2(press1)
	mod.temp <- lm(Chosen ~ backCapInt, weights=eWeight, data=temp)
	return(coef(mod.temp))
}

clusterExport(cl, list("press1", "clusterBootS2", "bootPress1", "bootPress1W"))

#Run models
system.time(boot.press1 <- parSapply(cl, rep(1,1500), bootPress1)) 
system.time(boot.press1W <- parSapply(cl, rep(1,1500), bootPress1W))

#Extract quantities of interest for main model
plot.press1 <- cbind(apply(boot.press1, 1, mean), apply(boot.press1, 1, quantile, c(0.025)), apply(boot.press1, 1, quantile, c(0.05)), apply(boot.press1, 1, quantile, c(0.95)), apply(boot.press1,1,quantile, c(0.975)))[-1,]

#Now add baseline categories for main model
plot.press1b <- rbind(rep(0,5), plot.press1)

#Extract quantities of interest for weighted model
plot.press1W <- cbind(apply(boot.press1W, 1, mean), apply(boot.press1W, 1, quantile, c(0.025)), apply(boot.press1W, 1, quantile, c(0.05)), apply(boot.press1W, 1, quantile, c(0.95)), apply(boot.press1W,1,quantile, c(0.975)))[-1,]

#Now add baseline categories for weighted model
plot.press1Wb <- rbind(rep(0,5), plot.press1W)

#Labels
pressLab <- c("Backed down, high capabilities, high stakes", "Backed down, low capabilities, high stakes", "Backed down, high capabilities, low stakes", "Backed down, low capabilities, low stakes", "Stood firm, high capabilities, high stakes", "Stood firm, low capabilities, high stakes", "Stood firm, high capabilities, low stakes", "Stood firm, low capabilities, low stakes")

##### Figure 3a ####
# Clean up labels and legend in post-production
dev.new(height=6.5, width=8.5)
par(mar=c(5.1,1.1,4.1,0), oma=c(0,4,0,0))
plot(plot.press1b[,1], nrow(plot.press1b):1, type="n", axes=FALSE, xlab="Average marginal component effect (AMCE)", ylab="",xlim=c(-0.6,0.2), ylim=c(0.75,8))

for (i in 1:nrow(plot.press1b)){
	points(plot.press1b[i,1], nrow(plot.press1b)-(i-1), pch=16)
	points(plot.press1Wb[i,1], nrow(plot.press1Wb)-(i-0.75), pch=12)
	lines(plot.press1b[i,c(2,5)], rep(nrow(plot.press1b)-(i-1),2))
	lines(plot.press1Wb[i,c(2,5)], rep(nrow(plot.press1Wb)-(i-0.75),2))
	lines(plot.press1b[i,c(3,4)], rep(nrow(plot.press1b)-(i-1),2), lwd=2)
	lines(plot.press1Wb[i,c(3,4)], rep(nrow(plot.press1Wb)-(i-0.75),2), lwd=2)				
	text(-0.3, nrow(plot.press1b)-(i-0.875), pressLab[i], cex=0.7)			
	}
abline(v=0)
axis(1)

### Alternative model specifications discussed in the caption of Fig 3a (plots not included in text)
press2 <- subset(press1, dat3$curBehavior=="Public threat") #Only cases w public threat
press3 <- subset(press1, dat3$curBehavior=="Public threat" | dat3$curBehavior == "Mobilized troops") #Only cases with public threat or mobilization

bootPress2 <- function(...){
	temp <- clusterBootS2(press2)
	mod.temp <- lm(Chosen ~ backCapInt, data=temp)
	return(coef(mod.temp))
}

bootPress3 <- function(...){
	temp <- clusterBootS2(press3)
	mod.temp <- lm(Chosen ~ backCapInt, data=temp)
	return(coef(mod.temp))
}

clusterExport(cl, list("press2", "press3", "bootPress2", "bootPress3"))

system.time(boot.press2 <- parSapply(cl, rep(1,1500), bootPress2))
system.time(boot.press3 <- parSapply(cl, rep(1,1500), bootPress3))

plot.press2 <- cbind(apply(boot.press2, 1, mean), apply(boot.press2, 1, quantile, c(0.025)), apply(boot.press2, 1, quantile, c(0.05)), apply(boot.press2, 1, quantile, c(0.95)), apply(boot.press2,1,quantile, c(0.975)))[-1,]

#Now add baseline categories
plot.press2b <- rbind(rep(0,5), plot.press2)

dev.new(height=6.5, width=8.5)
par(mar=c(5.1,1.1,4.1,0), oma=c(0,4,0,0))
plot(plot.press2b[,1], nrow(plot.press2b):1, type="n", axes=FALSE, xlab="Average marginal component effect (AMCE)", ylab="",xlim=c(-0.6,0.2))

for (i in 1:nrow(plot.press2b)){
	points(plot.press2b[i,1], nrow(plot.press2b)-(i-1), pch=16)
	lines(plot.press2b[i,c(2,5)], rep(nrow(plot.press2b)-(i-1),2))
	lines(plot.press2b[i,c(3,4)], rep(nrow(plot.press2b)-(i-1),2), lwd=2)		
	text(-0.3, nrow(plot.press2b)-(i-1), pressLab[i], cex=0.7)			
	}
abline(v=0)
axis(1)

plot.press3 <- cbind(apply(boot.press3, 1, mean), apply(boot.press3, 1, quantile, c(0.025)), apply(boot.press3, 1, quantile, c(0.05)), apply(boot.press3, 1, quantile, c(0.95)), apply(boot.press3,1,quantile, c(0.975)))[-1,]

#Now add baseline categories
plot.press3b <- rbind(rep(0,5), plot.press3)

dev.new(height=6.5, width=8.5)
par(mar=c(5.1,1.1,4.1,0), oma=c(0,4,0,0))
plot(plot.press3b[,1], nrow(plot.press3b):1, type="n", axes=FALSE, xlab="Average marginal component effect (AMCE)", ylab="",xlim=c(-0.6,0.2))

for (i in 1:nrow(plot.press3b)){
	points(plot.press3b[i,1], nrow(plot.press3b)-(i-1), pch=16)
	lines(plot.press3b[i,c(2,5)], rep(nrow(plot.press3b)-(i-1),2))
	lines(plot.press3b[i,c(3,4)], rep(nrow(plot.press3b)-(i-1),2), lwd=2)		
	text(-0.3, nrow(plot.press3b)-(i-1), pressLab[i], cex=0.7)			
	}
abline(v=0)
axis(1)

### Additional robustness check (footnote 22): does the interactive model fit the data better than an additive one?

mod.int <- lm(Chosen ~prevAct*Capabilities*Stakes, data=dat3)
mod.add <- lm(Chosen~ prevAct + Capabilities + Stakes, data=dat3)

anova(mod.int, mod.add) #p < 0.237 that difference in fit is significant
BIC(mod.int, mod.add)
22478.36-22511.59 #-33.23; additive model has better fit

##### Attribution theory hypothesis

mercer1 <- dat3[,-c(2:13,15:22,25)] #Slimmer dataframe

bootMercer1 <- function(...){
	temp <- clusterBootS2(mercer1)
	mod.temp <- lm(Chosen ~ desirable, data=temp)
	return(coef(mod.temp))
}

bootMercer1W <- function(...){
	temp <- clusterBootS2(mercer1)
	mod.temp <- lm(Chosen ~ desirable, weights=eWeight, data=temp)
	return(coef(mod.temp))
}

clusterExport(cl, list("mercer1", "bootMercer1", "bootMercer1W"))

system.time(boot.mercer1 <- parSapply(cl, rep(1,1500), bootMercer1))
system.time(boot.mercer1W <- parSapply(cl, rep(1,1500), bootMercer1W))

#Extract quantities of interest for main model
plot.mercer1 <- cbind(apply(boot.mercer1, 1, mean), apply(boot.mercer1, 1, quantile, c(0.025)), apply(boot.mercer1, 1, quantile, c(0.05)), apply(boot.mercer1, 1, quantile, c(0.95)), apply(boot.mercer1,1,quantile, c(0.975)))[-1,]

#Now add baseline categories
plot.mercer1b <- rbind(rep(0,5), plot.mercer1)

#Extract quantities of interest for weighted model
plot.mercer1W <- cbind(apply(boot.mercer1W, 1, mean), apply(boot.mercer1W, 1, quantile, c(0.025)), apply(boot.mercer1W, 1, quantile, c(0.05)), apply(boot.mercer1W, 1, quantile, c(0.95)), apply(boot.mercer1W,1,quantile, c(0.975)))[-1,]

#Now add baseline categories
plot.mercer1Wb <- rbind(rep(0,5), plot.mercer1W)

#Labels
mercerLab1 <- c("Ally, backed down", "Adversary, backed down", "Adversary, stood firm", "Ally, stood firm")

#Re-sort factors for easier visual comparison
mercerLab1 <- mercerLab1[c(1,4,2,3)]
plot.mercer1b <- plot.mercer1b[c(1,4,2,3),]
plot.mercer1Wb <- plot.mercer1Wb[c(1,4,2,3),]

#### Figure 3b #####
#Clean up labels and legend in post-production
dev.new(height=4, width=6.5)
par(mar=c(5.1,1.1,4.1,0), oma=c(0,4,0,0))
plot(plot.mercer1b[,1], nrow(plot.mercer1b):1, type="n", axes=FALSE, xlab="Average marginal component effect (AMCE)", ylab="",xlim=c(-0.3,0.2), ylim=c(0.75,4))

for (i in 1:nrow(plot.mercer1b)){
	points(plot.mercer1b[i,1], nrow(plot.mercer1b)-(i-1), pch=16)
	points(plot.mercer1Wb[i,1], nrow(plot.mercer1Wb)-(i-0.75), pch=12)
	lines(plot.mercer1b[i,c(2,5)], rep(nrow(plot.mercer1b)-(i-1),2))
	lines(plot.mercer1Wb[i,c(2,5)], rep(nrow(plot.mercer1Wb)-(i-0.75),2))
	lines(plot.mercer1b[i,c(3,4)], rep(nrow(plot.mercer1b)-(i-1),2), lwd=2)
	lines(plot.mercer1Wb[i,c(3,4)], rep(nrow(plot.mercer1Wb)-(i-0.75),2), lwd=2)				
	text(-0.3, nrow(plot.mercer1b)-(i-0.875), mercerLab1[i], cex=0.7)			
	}
abline(v=0)
axis(1)

### Alternative model specification discussed in the caption of Figure 3b (plot not included in text)
bootMercer2 <- function(...){
	temp <- clusterBootS2(mercer1)
	mod.temp <- lm(Chosen ~ desirable2, data=temp)
	return(coef(mod.temp))
}

clusterExport(cl, list("mercer1", "bootMercer2"))

system.time(boot.mercer2 <- parSapply(cl, rep(1,1500), bootMercer2)) 

#Extract quantities of interest
plot.mercer2 <- cbind(apply(boot.mercer2, 1, mean), apply(boot.mercer2, 1, quantile, c(0.025)), apply(boot.mercer2, 1, quantile, c(0.05)), apply(boot.mercer2, 1, quantile, c(0.95)), apply(boot.mercer2,1,quantile, c(0.975)))[-1,]

#Now add baseline categories
plot.mercer2b <- rbind(rep(0,5), plot.mercer2)

#Labels
mercerLab2 <- c("Ally, backed down against adversary", "Adversary, backed down against adversary", "Adversary, backed down against ally", "Adversary, stood firm against adversary", "Adversary, stood firm against ally", "Ally, backed down against ally", "Ally, stood firm against adversary", "Ally, stood firm against ally")

#Re-sort factors for easier visual comparison
mercerLab2 <- mercerLab2[c(1,6,7,8,2,3,4,5)]
plot.mercer2b <- plot.mercer2b[c(1,6,7,8,2,3,4,5),]

dev.new(height=6.5, width=6.5)
par(mar=c(5.1,1.1,4.1,0), oma=c(0,4,0,0))
plot(plot.mercer2b[,1], nrow(plot.mercer2b):1, type="n", axes=FALSE, xlab="Average marginal component effect (AMCE)", ylab="",xlim=c(-0.3,0.2))

for (i in 1:nrow(plot.mercer2b)){
	points(plot.mercer2b[i,1], nrow(plot.mercer2b)-(i-1), pch=16)
	lines(plot.mercer2b[i,c(2,5)], rep(nrow(plot.mercer2b)-(i-1),2))
	lines(plot.mercer2b[i,c(3,4)], rep(nrow(plot.mercer2b)-(i-1),2), lwd=2)		
	text(-0.3, nrow(plot.mercer2b)-(i-1), mercerLab2[i], cex=0.7)			
	}
abline(v=0)
axis(1)

#### Democratic credibility hypothesis

fearon <- subset(dat3, dat3$curBehavior!="Mobilized troops")[,-c(2:13,15:24)]

bootfearon1 <- function(...){
	temp <- clusterBootS2(fearon)
	mod.temp <- lm(Chosen ~ audCost, data=temp)
	return(coef(mod.temp))
}

bootfearon1W <- function(...){
	temp <- clusterBootS2(fearon)
	mod.temp <- lm(Chosen ~ audCost, weights=eWeight, data=temp)
	return(coef(mod.temp))
}

clusterExport(cl, list("fearon", "bootfearon1", "bootfearon1W"))

system.time(boot.fearon1 <- parSapply(cl, rep(1,1500), bootfearon1))
system.time(boot.fearon1W <- parSapply(cl, rep(1,1500), bootfearon1W))

#Extract quantities of interest for main model
plot.fearon1 <- cbind(apply(boot.fearon1, 1, mean), apply(boot.fearon1, 1, quantile, c(0.025)), apply(boot.fearon1, 1, quantile, c(0.05)), apply(boot.fearon1, 1, quantile, c(0.95)), apply(boot.fearon1,1,quantile, c(0.975)))[-1,]

#Now add baseline categories
plot.fearon1b <- rbind(rep(0,5), plot.fearon1)

#Extract quantities of interest for weighted model
plot.fearon1W <- cbind(apply(boot.fearon1W, 1, mean), apply(boot.fearon1W, 1, quantile, c(0.025)), apply(boot.fearon1W, 1, quantile, c(0.05)), apply(boot.fearon1W, 1, quantile, c(0.95)), apply(boot.fearon1W,1,quantile, c(0.975)))[-1,]

#Now add baseline categories
plot.fearon1Wb <- rbind(rep(0,5), plot.fearon1W)

fearonLab1 <- c("Dictatorship, nothing", "Democracy, nothing", "Democracy, public threat", "Dictatorship, public threat", "Mixed, nothing", "Mixed, public threat")

#Reorder factors for easier visual comparison
plot.fearon1b <- plot.fearon1b[c(1,4,5,6,2,3),]
plot.fearon1Wb <- plot.fearon1Wb[c(1,4,5,6,2,3),]
fearonLab1 <- fearonLab1[c(1,4,5,6,2,3)]

#### Figure 4 ####
# Clean up labels and legend in post-production
dev.new(height=5, width=6.5)
par(mar=c(5.1,1.1,4.1,0), oma=c(0,4,0,0))
plot(plot.fearon1b[,1], nrow(plot.fearon1b):1, type="n", axes=FALSE, xlab="Average marginal component effect (AMCE)", ylab="",xlim=c(-0.3,0.2), ylim=c(0.75,6))

for (i in 1:nrow(plot.fearon1b)){
	points(plot.fearon1b[i,1], nrow(plot.fearon1b)-(i-1), pch=16)
	points(plot.fearon1Wb[i,1], nrow(plot.fearon1Wb)-(i-0.75), pch=12)
	lines(plot.fearon1b[i,c(2,5)], rep(nrow(plot.fearon1b)-(i-1),2))
	lines(plot.fearon1Wb[i,c(2,5)], rep(nrow(plot.fearon1Wb)-(i-0.75),2))
	lines(plot.fearon1b[i,c(3,4)], rep(nrow(plot.fearon1b)-(i-1),2), lwd=2)	
	lines(plot.fearon1Wb[i,c(3,4)], rep(nrow(plot.fearon1Wb)-(i-0.75),2), lwd=2)			
	text(-0.25, nrow(plot.fearon1b)-(i-0.875), fearonLab1[i], cex=0.7)			
	}
abline(v=0)
axis(1)